Compatible Fossil Fuel CO 2 Emissions in the CMIP6 Earth System Models’ Historical and Shared Socioeconomic Pathway Experiments of the Twenty-First Century

: We present the compatible CO 2 emissions from fossil fuel (FF) burning and industry, calculated from the historical and Shared Socioeconomic Pathway (SSP) experiments of nine Earth system models (ESMs) participating in phase 6 of the Coupled Model Intercomparison Project (CMIP6). The multimodel mean FF emissions match the historical record well and are close to the data-based estimate of cumulative emissions (394 6 59 GtC vs 400 6 20 GtC, respectively). Only two models fall inside the observed uncertainty range; while two exceed the upper bound, ﬁve fall slightly below the lower bound, due primarily to the plateau in CO 2 concentration in the 1940s. The ESMs’ diagnosed FF emission rates are consistent with those generated by the integrated assessment models (IAMs) from which the SSPs’ CO 2 concentration pathwayswereconstructed;thesimplerIAMs’emissionsliewithintheESMs’spreadforsevenoftheeightSSPexperiments, theotherbeingonlymarginallylower,providingconﬁdenceintherelationshipbetweentheIAMs’FFemissionratesand concentration pathways. The ESMs require fossil fuel emissions to reduce to zero and subsequently become negative in SSP1-1.9, SSP1-2.6, SSP4-3.4, and SSP5-3.4over. We also present the ocean and land carbon cycle responses of the ESMs in the historical and SSP scenarios. The models’ocean carbon cycle responsesare in close agreement, but there is considerable spread in their land carbon cycle responses. Land-use and land-cover change emissions have a strong inﬂuence over the magnitude of diagnosed fossil fuel emissions, with the suggestion of an inverse relationship between the two.


Introduction
The release of carbon dioxide and other greenhouse gases since the onset of the industrial era has warmed the planet by around 1 18C (Stocker et al. 2013).The changing climate is already affecting natural terrestrial and marine ecosystems, thawing permafrost, and altering the hydrological cycle, with implications for drought, flooding, and wildfires.Human systems are being impacted increasingly, with reductions in crop yields, disruption to infrastructure, water supply and food production, and increased mortality due to high and low temperature extremes (IPCC 2014).With global emissions of CO 2 continuing to rise (Friedlingstein et al. 2019), the 2015 Paris Agreement of the United Nations Framework Convention on Climate Change (UNFCCC) saw 196 governments around the world agree to a commitment to constrain greenhouse gas (GHG) emissions to a level that would limit long-term warming, relative to the preindustrial level, to well below 28C, ideally to below 1.58C, with emissions peaking as soon as possible and then reducing to zero within decades (UNFCCC 2015).
In its Fifth Assessment Report, the Intergovernmental Panel on Climate Change (Collins et al. 2013) reported for the first time that the change in global mean temperature relative to the preindustrial is proportional to the cumulative total CO 2 emissions released, albeit with significant uncertainty.Whether warming is limited to 1.58 or 28C above preindustrial levels is therefore greatly dependent on how CO 2 emissions proceed over the coming century.How society evolves in the future, how the global population will change, and the extent to which countries and communities will cooperate to reduce emissions from the energy sector, transport, and agriculture cannot be known with any certainty.However, integrated assessment models (IAMs), computer models of intermediate complexity, can be used to estimate how societal changes such as population growth, economic development, investment priorities, and degree of international cooperation will impact future greenhouse gas emission rates and concentrations (Gidden et al. 2019).
Phase 6 of the Coupled Model Intercomparison Project (CMIP6) is divided into several individual model intercomparison projects (MIPs), each with a particular science focus (Eyring et al. 2016), undertaking experiments with fullcomplexity Earth system models (ESMs) and general climate models.ScenarioMIP (O'Neill et al. 2016) is designed to compare ESMs' responses to eight new forcing scenarios out to 2100 built around five Shared Socioeconomic Pathway (SSP) scenarios of potential future social and economic development (Riahi et al. 2017).
The current generation of Earth system models typically contain fully process-based representations of the terrestrial and ocean carbon cycles, allowing uptake of carbon by land and the oceans in response to changes in atmospheric CO 2 concentration to be quantified.These fluxes can be used in conjunction with the prescribed changes in atmospheric CO 2 burden to determine the rate of fossil fuel emissions that are compatible with the prescribed CO 2 concentration pathway (see section 2a).When an ESM is forced with observed historical CO 2 concentration, the compatible fossil fuel (FF) emissions can be compared with historical anthropogenic FF emissions as a test of the model's ability to reproduce various components of the historical carbon budget realistically.For future scenarios, compatible emissions from ESMs can be used to evaluate the much less complex integrated assessment models; a good correlation between IAM emissions and ESM emissions provides evidence that the IAM's carbon cycleclimate submodel, which generates the CO 2 concentration pathway from its emissions, is broadly consistent with the ESM's much more sophisticated carbon cycle.While an IAM provides one ''realization'' of future emissions compatible with a CO 2 pathway, with the different responses of their land and ocean carbon cycle components, ESMs provide a range of emissions that are compatible with a given CO 2 pathway.This range highlights the uncertainty in our understanding of the response of the coupled climate-carbon cycle system.It provides, therefore, an estimate of uncertainty in the compatible emissions rate, cumulative emissions, and the FF contribution to associated carbon budgets; for example, emissions remaining before 1.58 or 2.08C of warming are reached.
In section 2, we discuss in more detail the way in which the compatible emissions are diagnosed, as well as the SSP scenarios.We present the compatible emissions from the nine ESMs' historical and SSP experiments in section 3 and examine which scenarios require emissions to reduce to below zero, and the timing at which this occurs.We show the ESMs' land and ocean carbon cycle responses, as well as the cumulative airborne fraction of emissions.The nine CMIP6 ESMs included in this study are summarized in Tables 1-3 and are described in further detail in section 1 of the online supplemental material.

a. Calculation of compatible fossil emissions
Although the current generation of Earth system models typically contain a fully interactive carbon cycle in which the model can simulate the atmospheric CO 2 concentration in response to a scenario of CO 2 emissions, when taking part in model intercomparison projects, such models typically use prescribed CO 2 .This makes it simpler to ensure that the resulting spread of climate projections are in response to identical CO 2 forcing.However, since the atmospheric CO 2 is prescribed, and the land and ocean carbon cycle submodels calculate how much CO 2 is taken up by land and by the oceans, it is possible to infer the rate at which additional CO 2 must be emitted to maintain the evolving atmospheric CO 2 pathway, offsetting that taken up by the land and ocean sinks.
Implicit N limitation: LAI reduces as CO 2 increases, along with downregulation of photosynthesis As described in appendix 1 of Jones et al. (2013), any increase to the combined atmosphere, land and ocean carbon store over time must be equal to the flux of ''new'' carbon added to the system by the combustion of fossil fuels and other processes that release to the atmosphere carbon that was previously locked away on geological time scales, such as gas flaring and cement manufacture (when using the term ''fossil fuel,'' we include all of these sources).By contrast, anthropogenic emissions from land-use practices release to the atmosphere CO 2 that is continually cycling through the land, ocean, and atmosphere carbon stores, on time scales from weeks to centuries, and so represent no net increase to the combined global carbon stores.
Equation (1), below, adapted from Eq. ( 1) of Friedlingstein et al. (2019) shows the relation between growth of the atmospheric carbon store (G ATM ), the gross land carbon sink before accounting for emissions due to land-use change (S LAND ), the ocean carbon sink (S OCEAN ), and anthropogenic emissions of CO 2 from fossil sources (E FF ), and from land-use change (E LUC ): (1) The budget imbalance term of Eq. ( 1) of Friedlingstein et al. (2019) is a residual term required to close the budget since all other terms are based on observations, and are therefore inexact.In the context of an individual model's simulated carbon stores and the fluxes between them, all terms are known exactly, so an imbalance term is not necessary.Rearranging Eq. ( 1) gives We calculate the annual atmospheric CO 2 growth rate (G ATM ) in gigatons of carbon (1 GtC 5 1 PgC 5 10 15 gC) per year, from the prescribed CO 2 concentration in parts per million (ppm) using the equivalence 1 part per million of CO 2 in the atmosphere is equal to 2.124 GtC (Ballantyne et al. 2015).The annual mean ocean carbon sink (S OCEAN ) is calculated from the global total of the air-to-sea flux of CO 2 .The annual mean net land sink after accounting for emissions from land use (S LAND 2 E LUC ) is equal to the global total of the net biosphere production: the gross land carbon sink (S LAND ) reduced by emissions due to land-use change, harvest, grazing, and fire (E LUC ).The net biosphere production (NBP) diagnostic, nbp, and the air-to-sea CO 2 flux diagnostic, fgco2, are two of the variables submitted to the Earth System Grid Federation (ESGF) CMIP6 data archive.2Thus, substituting the annual mean growth rate of the prescribed CO 2 burden, fgco2, and nbp into (2) we calculate the annual mean fossil fuel emissions rate (E FF ) shown in section 3.

b. Calculation of emissions due to land-use and land-cover change
Emissions of CO 2 associated with wood harvesting and landuse and land-cover changes (LULCC) represent the greatest source of uncertainty in the global carbon cycle (Friedlingstein et al. 2019).These emissions cannot be measured directly; observations can be made only of the net flux of CO 2 between the land and the atmosphere, that is, the natural land carbon sink minus carbon lost due to LULCC emissions.Isolating the latter and its associated uncertainty has been the focus of enormous effort for many years (Pongratz et al. 2014).The last two Global Carbon Budget estimates of LULCC emissions (Le Quéré et al. 2018;Friedlingstein et al. 2019) have been calculated as the mean of the results from two bookkeeping models (Houghton and Nassikas 2017;Hansis et al. 2015).Bookkeeping models first consider literature-derived estimates of the changes in soil and vegetation carbon stores, before and after a change in land use such as replacement of forest with crops.They then distribute the resulting carbon stock changes in time according to observation-based curves of decay of plant and soil carbon in the field or in products, and regrowth following agricultural abandonment or wood harvesting.
A multitude of terms has been used in the literature to refer to emissions from land use, with contributions from a great range of constituent processes.The scope of processes that can be simulated depends on the nature of the type of model in question, with bookkeeping models, dynamic vegetation models and fully coupled Earth system models differing in their capacity to include them (Pongratz et al. 2014).In the experiments considered here, since the atmospheric CO 2 concentration is prescribed, the CO 2 emitted by LULCC is not added to the atmosphere store.
The nine ESMs considered in the present study all include some representation of land-use change, of varying degrees of complexity, with associated CO 2 emissions.Land management processes, such as crop harvest, wood harvest, and grazing of pasture, are included in some models but not others, as summarized in Tables 1-3.The land-use emissions in a single experiment will include the contribution from the processes specific to each model, output by the appropriate combination of diagnostics; this approach to quantifying LULCC emissions is described as method S (single simulation) in Pongratz et al. (2014).
A more accurate method of determining the LULCC emissions E LUC is to calculate the difference in net biosphere production between a pair of simulations, one with land use changing over time, and the other with fixed land use [the experimental design denoted E2 in Pongratz et al. (2014)].For the historical period, this is given in Eq. ( 3), where the two experiments are the standard CMIP6 historical, and the fixed, preindustrial land-use variant, ''hist-noLu'', included in the Land Use Model Intercomparison Project (LUMIP; Lawrence et al. 2016): (3) This approach identifies the total difference in land carbon store due solely to land-use change and associated processes in the historical experiment.It takes into account the additional carbon uptake as natural vegetation is allowed to regrow on abandoned agricultural land, which is neglected by the online diagnosis of land-use emissions within a single simulation.It also accounts for the loss of the additional sink capacity (LASC) that natural vegetation would have provided had it not been replaced by managed land.The inclusion of the loss of additional sink capacity constitutes a key difference to estimates from the bookkeeping approach, which does not consider transient changes in environmental conditions (Pongratz et al. 2014).By using only the diagnostic of net biosphere production nbp, the ''E2'' approach does not rely on identifying the combination of diagnostics specific to each model that constitute its total LULCC emissions.However, it does require the additional computational expense of performing the no-land-use change noLu variant of experiments with landuse change.

c. Scenarios
Under the CMIP5 modeling framework (Taylor et al. 2012), the response of the Earth system to realistic scenarios of future forcing was explored via the representative concentration pathways (RCPs; Vuuren et al. 2011).These consist of four distinct pathways leading to radiative forcing at 2100 of approximately 2.6, 4.5, 6.0, and 8.5 W m 22 , known as RCP2.6,RCP4.5, RCP6.0, and RCP8.5, respectively.Each pathway provides the greenhouse gas concentration and land-use scenarios, such as the area of land devoted to crop and pasture, required as forcing by Earth system models.
For CMIP6, there are eight future scenario experiments that make up ScenarioMIP (O'Neill et al. 2016); four share similar 2100 forcing levels with the RCPs, while four more lead to complementary radiative forcing levels of 1.9, 3.4 (two scenarios), and 7.0 W m 22 .The basis for each scenario is one of five SSPs, numbered 1 to 5 (Riahi et al. 2017).The SSPs describe five future narratives of economic development (Dellink et al. 2017), population change (Samir and Lutz 2017), and urbanization (Jiang and O'Neill 2017).They are designed to encompass a range of potential futures characterized by how challenging the mitigation of climate change, and adaptation to it, will be in those futures.
The scenarios are generated by intermediate-complexity integrated assessment models, which often comprise many individual submodel components (Riahi et al. 2017).They take as their inputs a wide range of variables such as scenarios of economic development, population growth, the degree of international cooperation, policy decisions on the importance of investing in education, reducing inequality, sustainability of energy resources and many other complex traits of modern society.The models are run into the future to determine how energy use and emissions of CO 2 and other greenhouse gases evolve in response to these input scenarios.For ScenarioMIP, the emissions of different greenhouse gases from the IAMs were converted into their respective concentrations using version 7.0 of MAGICC, a carbon cycle-climate model of intermediate complexity (Meinshausen et al. 2020(Meinshausen et al. , 2011)).
The ScenarioMIP experiments are labeled SSPx-y.zwhere x is the SSP number from 1 to 5, and y.z is the forcing level in W m 22 , as was the case with the RCPs.The CMIP6 experiment identifier for each of the SSP experiments is sspxyz, for example, ssp370, with the exception of ssp534-over, which is an overshoot scenario, branching from ssp585 in 2040.The experiments are ordered by priority into tier 1 (ssp126, ssp245, ssp370, ssp585) and tier 2 (ssp119, ssp434, ssp460, ssp534-over).
We present results from all eight SSP experiments, whose CO 2 concentration pathways are shown in Fig. 1a, with the global population in Fig. 1b.The crop and pasture fractions of the land-use scenarios are shown in Fig. 1c, with forest area in Fig. 1d.

e. Reference data
Where possible, the historical plots and tables include reference data from the Global Carbon Budget, 2019 (GCB2019; Friedlingstein et al. 2019).The global fossil emissions rate (from fossil carbonates as well as fossil fuels) E FF are from Gilfillan et al. (2019).The E LUC emissions from GCB2019 are the average of updated versions of the data from Houghton and Nassikas (2017) and Hansis et al. (2015).The terrestrial and ocean carbon sinks, S LAND and S OCEAN , respectively, come directly from the GCB2019.The E FF is based on historical records of fossil emissions, but the E LUC , S OCEAN , and S LAND are all derived from models constrained by observations where possible.

Results and discussion
The results have been divided into a historical section from 1850 to 2014 (section 3a) and a future section covering all SSP experiments from 2015 to 2100 (section 3b).The time series shown were constructed from ensemble means, where possible, of the nbp, fgco2, cVeg, cSoil, cLitter (vegetation, soil, and litter carbon store), fharvest (harvest flux), and tas (surface air temperature) diagnostics submitted to ESGF.The models' individual responses to the SSPs for each variable are shown in the online supplemental material, which also contains tables summarizing ocean and land carbon uptake historically and under the SSPs.

1) COMPATIBLE FOSSIL FUEL EMISSIONS
Equation (2) has been used to calculate time series of compatible fossil fuel CO 2 emissions for the ESMs' historical experiment (Fig. 2).For comparison, the data-based estimates of historical emissions from fossil fuel burning and cement production (Gilfillan et al. 2019) are shown (black dashed lines).The multimodel mean follows very closely the observed historical emission rate, except for a period of around 15 years, from the early 1940s.After decades of rising, in 1941 CO 2 began to stabilize and remained between 310 and 312 ppm from 1941 to 1951 before rising again [Etheridge et al. (1996), confirmed more recently by MacFarling Meure et al. (2006) and Rubino et al. (2013)].Hence, during this period of the historical experiment, G ATM reduced to almost zero, leading the ESMs to simulate a reduction in E FF , in the absence of any compensating factor to increase the ocean or land sinks.In reality, the data-based fossil fuel emissions rate showed no such lull during the observed plateau in concentration at this time, increasing from 1.33 to 1.76 GtC yr 21 from 1941 to 1951.This implies the existence in the real world of an enhanced carbon sink during this period that has not yet been explained fully by the observational datasets available (Bastos et al. 2016).The ocean is likely to have played a role, with enhanced uptake by the Southern and tropical Pacific Oceans associated with the strong El Niño of the early 1940s (Joos et al. 1999).There could have been additional uptake on land due to land-use processes and changes not accounted for in the LUC datasets, nor therefore in estimates of land uptake of the time.Due to the 1-2.6 2-4.5 3-7.0 4-3.4 4-6.0 5-3.4over 5-8.5 Second World War, the 1940s were characterized by enormous socioeconomic upheaval with changes to land management, food production, and relocation of population; not all changes contributing to LUC emissions were recorded accurately, such as abandonment of land and changes to wood harvest (Bastos et al. 2016).It should also be noted that this era predates direct measurement of atmospheric CO 2 concentration, and estimates of fossil fuel emissions were less accurate than they are today.So our understanding of the 1940s CO 2 plateau and of the fluxes potentially contributing to the implied sink remains incomplete (Bastos et al. 2016).It would be wrong to conclude, therefore, that because the ESM's fossil emissions fall short of the observations-based emissions during this period, the models systematically underestimate FF emissions during periods of constant or very slowly changing CO 2 , as occur in some of the SSP scenarios discussed in section 3b.
The ESMs' temperature response to the historical forcing is shown in Fig. S1 in the online supplemental material.

2) CUMULATIVE FOSSIL FUEL EMISSIONS
Figure 2b shows the cumulative fossil fuel emissions from the ESMs, with the data-based estimate shown for comparison (Gilfillan et al. 2019).The mean and standard deviation of the total FF emissions from 1850 to 2014 was 394 6 59 GtC, compared to observations-based estimate of 400 6 20 GtC for the same period (Gilfillan et al. 2019).The model mean is therefore very close to the observed estimate, lying well within its uncertainty range (Table 5); however, exclusion of the highest two model values results in a model mean of 367 GtC, below the lower limit (380 GtC) of the observed range.The ESMs' underestimate of emissions during the period of CO 2 plateau is largely responsible for this; cumulative totals of  several of the ESMs diverge from the observations-based total at this time and remain so until 2014 (Fig. 2b).Only two of the nine ESMs (IPSL and MIROC) lie within the observed cumulative 1850-2014 range (Table 5).Before looking at how the land and ocean sinks contribute to the historical FF emission rate for the nine ESMs in detail, it is worth examining a subset of simulations with fixed land use, to see the extent to which land-use change influences the magnitude of, and model spread in, compatible fossil fuel emissions.

3) INFLUENCE OF HISTORICAL LAND-USE EMISSIONS ON FOSSIL FUEL EMISSIONS
Of the nine ESMs included in this study, seven have submitted output from both historical and hist-noLu experiments to the CMIP6 archive.The LULCC emissions rate, calculated using Eq. ( 3), is shown in Fig. 3c, with the cumulative total in Fig. 3d, for comparison with the equivalent FF emissions (Figs.3a,b, respectively) from the historical, with GCB2019 reference values for comparison.Note that the GCB2019 LULCC emission estimate is based on the bookkeeping approach, which does not include loss of additional sink capacity.The LASC adds approximately another 20% to the LULCC emissions estimate [GCB2019 states 10%, but this was revised in Gasser et al. (2020)].For clarity, larger versions of the individual panels of Fig. 3 appear in the supplemental material as Figs.S13-S15.
The spread in the magnitude of the LULCC emissions among the seven models is immediately apparent, with CESM2's, which are very close to the observed estimates, being more than eight times those from CNRM.However, it is noticeable that all ESMs exhibit lower cumulative E LUC than GCB2019, with the exception of the last few years of MPI, despite their accounting for the loss of additional sink capacity.This hints toward LUC processes missing from the ESMs that are included in the GCB2019 estimate due to the construction of the bookkeeping models, which use observed carbon densities for various managed and unmanaged vegetation types.
The model spread in LULCC emissions is mirrored approximately by the spread in FF emissions, cumulatively (Figs. 3d,b,respectively); greater LULCC emissions are paired with lesser FF emissions, as is seen most clearly with those ESMs at the edges of the model spread, CNRM, MPI, and CESM2.Broadly speaking, the greater the extent to which a model's land carbon uptake is limited by land use, the lesser the amount of fossil CO 2 can be emitted for a given CO 2 concentration pathway.Combining the FF and LULCC emissions (Fig. 3e) significantly reduces the model spread (Fig. 3f) as the influence of land use on both LULCC and FF emissions cancels out, so differences only in the models' natural carbon sink are in evidence; the model spread at 2013 for FF 1 LULCC is 77 GtC, compared to 168 GtC for FF and 173 GtC for LULCC.The differences in the models' natural carbon sink are therefore dwarfed by the differences between the net land carbon sink after accounting for land-use practices (as will be seen in Fig. 5b).For the seven models included in Fig. 3, the presence of land use is therefore the greatest source of uncertainty in future FF emissions, and on future carbon budgets utilizing these.
Similar analysis cannot be done for the SSP experiments since the noLu variant of the experiments have not been requested under CMIP6 protocols for ScenarioMIP, C4MIP or LUMIP.
The magnitude of LULCC emissions calculated from the E2 paired experiments depends on several factors, including the spatial extent and dynamics of the land-use forcing dataset used in the historical; the carbon density of the vegetation being cleared; the land-use practices such as irrigation, fertilization, harvesting, and grazing, which impact the agricultural vegetation; and the responses of both the natural and agricultural vegetation to the changing climate and CO 2 concentration.The response of natural vegetation to future climate and CO 2 will depend on whether the vegetation distribution is fixed, or a dynamic vegetation model is utilized that allows the vegetation cover to migrate to latitudes and regions that become habitable as the climate changes.Land management practices that transfer living biomass to product pools, such as wood harvest and crop harvest, also play a role.For the historical experiment, although the spatial land-use forcing is provided in each case by the Land Use Harmonization 2 (LUH2) dataset (Hurtt et al. 2020), the combination of landuse states from LUH2 used is model dependent.So too is the method by which vegetation types within a grid cell change.''Net transitions,'' considering only the net differences of each vegetation type's fractional coverage following a land-cover change, can underestimate the impact on carbon stores compared to the ''gross transitions'' approach.Under gross transitions, the localized carbon impact at subgrid scale of each individual transition between neighboring vegetation types is considered explicitly; for example, primary forest of a young age structure replaces abandoned cropland, which is distinct from older, secondary forest already present elsewhere in the gridbox (Yue et al. 2018).Of the nine ESMs in the present study, two (MIROC, MPI) consider gross transitions; these, along with CESM2, get closest to the GCB2019 cumulative E LUC .

4) LAND CARBON RESPONSE
The simulated historical nbp, the land carbon uptake after accounting for carbon lost due to land-use change, is shown in Fig. 4. The GCB2019 provides observationally constrained multimodel estimates of the land carbon sink in each of the decades from the 1960s onward; these are included in Fig. 4 and are compared to the model's means for the same time periods in Table S2.The model mean sits well within the uncertainty range of the decadal mean estimates.The GCB2019s S LAND and E LUC estimates of 185 6 50 GtC and 195 6 75 GtC, respectively, result in net land carbon uptake from 1850 to 2014 of 210 6 125 GtC.The ESMs' net uptake ranges from 235 to 1142 GtC over this period, with a mean of 120 GtC, as summarized in Table S3.
The nine ESMs' preindustrial (1850) land carbon stores, with partitioning between soil, litter, and vegetation, are shown in Fig. 5a.The change in land carbon (cumulative nbp) relative to 1850 is shown in Fig. 5b, with changes in vegetation carbon (cVeg) and soil carbon (cSoil 1 cLitter) in Figs.5c and 5d, respectively.Broadly, the models' historical land carbon uptake is dominated by the change in vegetation carbon, with the initial reduction due to deforestation offset in later decades by The multimodel mean change in soil carbon shows a very gradual increase throughout the historical period.By contrast, MIROC's soil carbon reduces until the final few decades, when it begins to increase.In MIROC, replacement of primary/ secondary vegetation by crop and pasture of lower productivity reduces the soil and litter carbon until they are balanced by the lower litter inputs.In the latter half of the twentieth century, due to the effects of CO 2 fertilization and nitrogen fertilizer application, vegetation, litter, and soil carbon increase.Considering MIROC's intermediate CO 2 concentration-carbon feedback strength (Arora et al. 2020), the relatively steep increase of land carbon in this period is likely induced by accumulation in cropland, as suggested in Hajima et al. (2020).
The presence of a crop harvest flux can impact the vegetation carbon and/or soil carbon stores significantly.The crop harvest from five ESMs' historical simulations relative to their piControl is shown in Fig. 6, as a rate and cumulatively.The cumulative total at 2014 ranges from 60 to 246 GtC, so represents a fairly large impact on vegetation, litter, or soil carbon stores depending on how harvest is configured in each model.In UKESM, for example, 30% of the litter flux from crops to soil is intercepted to represent harvest, so UKESM's increase in historical soil carbon is not as great as it would otherwise be, but its living vegetation carbon store is not directly affected by the harvest.
Of the three ESMs that do not include a crop harvest, two (ACCESS, CNRM) exhibit the two greatest increases in historical soil and land carbon, contributing to the two greatest fossil fuel emission rates of all nine ESMs.Conversely, the third ESM in which crop harvest is not represented, CanESM5, sees a downward trend in soil carbon throughout the historical period that steepens in the final few decades.This is because in CanESM5 the soil decomposition rate over croplands is higher and the fraction of humified litter that is transferred to the soil carbon, as opposed to respired to the atmosphere, is lower than over natural vegetation.Consequently, as natural vegetation is replaced by croplands over the historical period, a decrease in global soil carbon is obtained, as is also seen in empirical measurements (Wei et al. 2014).
During harvest in CESM2 and NorESM2 (in both of which the land component is CLM5), even though crop grain carbon is removed, the increased productivity of crops due to fertilization, irrigation, and higher productivity plants actually leads to increases in soil carbon where crops are grown.This response is the opposite to what is observed in the real world, which is possibly due to lack of a representation of tillage in the model; tillage would likely lead to increased respiration, resulting in soil carbon losses.Further, what happens with soil carbon will depend on how much plant matter is removed during harvest; CLM5 assumes only grain is removed, but in reality, usually more than just grain is removed during the harvest process (Lombardozzi et al. 2020).
CESM2 and NorESM2 also simulate a wood harvest flux, which can result in a greater loss of vegetation carbon than deforestation.The wood harvest flux has a greater effect on vegetation carbon store as the regrowth of the harvested trees occurs on a much longer time scale than is the case with crops.This also impacts the soil carbon evolution in CESM2.
MIROC and MPI include an additional carbon loss due to grazing of pasture, which in MIROC results in a reduction in leaf area and reduces productivity.

5) OCEAN CARBON RESPONSE
The air-to-sea CO 2 fluxes are very similar in all nine ESMs (Fig. 7a) throughout the historical period.The model mean lies slightly lower than the six decadal mean values from GCB2019, but well within their uncertainty range, as is true for all individual models; these are compared with the GCB2019 reference values in Table S4.The models sit at (ACCESS) or below the GCB2019 change in ocean carbon store of 150 GtC from 1850 to 2014 (Table S4).IPSL and CNRM fall just below the lower end of the uncertainty range of the reference values from GCB2019.

6) AIRBORNE FRACTION
The rate of growth of atmospheric CO 2 concentration is dependent on the airborne fraction (AF) of emitted CO 2 , the proportion of emissions that remains in the atmosphere rather than being absorbed by the oceans or land.The AF varies with the rate of emissions; if uptake by the sinks is outpaced by emissions, the emitted CO 2 will begin to accumulate in the atmosphere, resulting in a higher AF and greater radiative forcing.Whether the airborne fraction of emissions has begun  to increase in recent history and whether it will do so over the coming century are therefore important questions.
Several studies in recent years have attempted to find evidence for a trend in AF from observations taken over the last 50 years or so.The airborne fraction is typically determined as a fraction of total anthropogenic emissions, from fossil fuel as well as land-use change (Canadell et al. 2007;Denman et al. 2007;Le Quéré et al. 2009;Raupach et al. 2008;Rayner et al. 2010;Bennedsen et al. 2019), though it has also been estimated from fossil fuels and other industrial sources alone, for example, from observations by Keeling et al. (1995) and for the CMIP5 models by Jones et al. (2013) and Ciais et al. (2013).These studies have suggested a weak increase in the airborne fraction of emitted carbon since the 1960s, although the robustness of this finding has been questioned (Knorr 2009;Raupach et al. 2014;Ballantyne et al. 2015).Since the focus of this paper is on compatible fossil fuel emissions, and due to the range of approaches to land-use emissions in the CMIP6 ESMs, we calculate the AF from the fossil fuel emissions alone, and the cumulative airborne fraction (CAF) rather than the instantaneous annual mean AF.The annual mean AF is a quotient constructed from three potentially small fluxes (rise in atmospheric CO 2 , nbp, and fgco2 for a particular year) and can vary widely from year to year, whereas the cumulative AF is calculated as the total amount by which the atmospheric burden has increased since the preindustrial, expressed as a fraction of the cumulative fossil fuel emissions.The cumulative AF has a much stronger signal-to-noise ratio after several decades of emissions and will change over time to reflect the ability of the land and ocean sink to keep pace with the emissions.
The time series of cumulative airborne fraction for the ESMs in the historical period is shown in Fig. 8.There is little  agreement between models in the trend and magnitude of the CAF historically.Four models' FF emissions are exceeded by the increment in atmospheric CO 2 resulting in a CAF exceeding unity, until the final few decades of the twentieth century; as a result, the same is true of the model mean CAF until 1940.By the end of the twentieth century, all models are converging toward a CAF of between 0.6 and 0.7, with the exception of the two ESMs in which the land carbon store increases throughout the historical period (ACCESS and CNRM).The IPCC Fifth Assessment Report gives estimates of the airborne fraction of fossil fuel emissions for the 1980s, 1990s, and two decadal periods of the present century, which also exhibit a reducing trend, from 0.62 in the 1980s to 0.51 in the first decade of this century (black boxes with uncertainty range on Fig. 8).

7) POTENTIAL ROLE OF NON-CO 2 GHGS AND AEROSOLS ON LAND AND OCEAN CARBON UPTAKE
Carbon uptake by the biosphere can be influenced directly by forcings other than CO 2 (Gasser et al. 2017).Aerosolinduced cooling can impact land carbon uptake directly.Zhang et al. (2019) showed that aerosol cooling increased historical land carbon uptake by 11.6 to 41.8 GtC, due primarily to increases in net biosphere production in the mid-to high latitudes; at lower latitudes, aerosol cooling reduced gross primary production (GPP) by more than total ecosystem respiration, reducing NBP.
Complex aerosol-cloud interactions are expected to lead to changes in precipitation, with implications for land carbon uptake (e.g., Fig. S11 of Zhang et al. 2019).Unger et al. (2017) used a coupled aerosol-carbon cycle-climate model to simulate the period 1996 to 2005 with and without historical aerosol levels and found that the inclusion of aerosols led to a global precipitation reduction of 0.08 mm day 21 with regional increases and decreases in NPP.
A further influence of aerosols on land carbon uptake is via their increased scattering of light, increasing the proportion of diffuse to direct light and allowing it to penetrate further into the canopy, enhancing photosynthetic capacity and increasing GPP (Mercado et al. 2009).
ESMs differ markedly in their representation of aerosols and atmospheric chemistry of aerosol precursor species (e.g., Thornhill et al. 2021).It is beyond the scope of this paper to examine in depth these differences and the resulting aerosol concentrations, but the effective radiative forcing (ERF) 3 from 1850 to 2014 due to the presence of aerosols (ERF aer ) has been calculated by Smith et al. (2020) for 6 ESMs in the present study.The ERF aer ranges from 20.63 W m 22 (IPSL), to 21.37 W m 22 (CESM2), a span of 0.74 W m 22 .This illustrates that aerosol-induced climate cooling and associated changes to precipitation could vary quite strongly between models, with the land carbon uptake, and therefore the compatible fossil emissions, influenced accordingly.
The AerChemMIP (Collins et al. 2017) experimental design allows calculation of the ERF at 2014 due to well-mixed greenhouse gases (CO 2 , CH 4 , N 2 O, and halocarbons, ERF wmghg ), CO 2 only (ERF CO2 ), and, therefore, well-mixed GHGs excluding CO 2 (ERF nonCO2ghg ).The range of ERF nonCO2ghg is smaller than ERF aer , from 0.71 to 1.15 W m 22 , a span of 0.44 W m 22 ; this is smaller than the spread in ERF aer but shows that differences in the models' responses to the non-CO 2 greenhouse gas forcing also plays a nontrivial role in the evolution of the ESMs' climate and therefore of the carbon cycle and compatible emissions.Ozone can also affect land carbon sinks directly, through its impact on leaf health (Sitch et al. 2007).
Aerosol-induced cooling also increases solubility of gaseous CO 2 in seawater, and to some extent alters ocean circulation and biological production, all of which potentially impact ocean carbon sinks, albeit to a lesser extent than that expected 3 ERF is the radiative forcing (RF) after allowing for the adjustment of all tropospheric and land surface properties on a range of time scales (though typically referred to as ''rapid adjustments''), including changes to vegetation and land ice cover, and is evaluated after these have regained equilibrium following a perturbation in forcing such as a change in CO 2 or aerosol.By contrast, RF is defined as the change in downward radiative flux at the tropopause after allowing only stratospheric temperatures to adjust to radiative equilibrium, with all other state variables held constant at the unperturbed state.The ERF represents a better indicator than RF of the eventual global mean temperature in response to a perturbation in forcing.For aerosols, in particular, ERF can be quite different from RF (Smith et al. 2020).
Unauthenticated | Downloaded 05/20/24 03:33 AM UTC from the land biosphere (Tjiputra et al. 2016;Lauvset et al. 2017).Several ESMs (MIROC, NorESM2, and UKESM) further simulate oceanic dimethyl sulfide (DMS) emissions, which is the largest natural source of atmospheric sulfur and contributes to aerosol formation and can therefore feed back to the oceanic carbon sink.
b. Future

1) COMPATIBLE FOSSIL FUEL EMISSIONS
For the SSP experiments, the ESM FF emissions are compared with the emission rates from the IAMs in Fig. 9a.There is generally strong agreement between the IAM and the ESMs; the IAM FF emission rates lie within the spread of the ESMs, with the exception of SSP4-6.0, which falls slightly below the ESMs' range.The SSP4-6.0 was performed by fewer ESMs (three) than any other SSP experiment; these include CNRM and CanESM, which simulate two of the three largest cumulative FF emissions for SSP2-4.5, the scenario with emissions closest to those of SSP4-6.0.Had more ESMs performed SSP4-6.0, the range of ESM emissions might well have encompassed those of the IAM.(The fossil fuel emission rates are shown for the models individually in Fig. S3.) Three of the SSPs require the emissions to reduce to below zero by 2100 to adhere to the CO 2 pathway, in at least one ESM.In all ESMs, SSP1-2.6 emissions become negative, though only very slightly in NorESM2, and the same is true of the five models that performed SSP1-1.9 and SSP4-3.4.The interannual variability of compatible emissions varies a lot between ESMs; some have very little (CanESM5, CNRM, CESM2), others more so (ACCESS, MIROC).MIROC's emissions are characterized by a fairly cyclical variability with a period of about 10 years and an amplitude of 5 GtC yr 21 .This is likely due to the unrealistically large El Niño-Southern Oscillation amplitude exhibited by MIROC that causes similar variability in global temperature (Hajima et al. 2020) as can be seen in Fig. S2.
Of interest to policymakers, particularly in the case of overshoot scenarios, is the point at which emission rates peak and start to decline.Figure 10 shows the peak emissions rate versus the year in which this occurs for all nine ESMs and the IAMs.For the lower emissions scenarios, the IAM and ESMs are in close agreement for the magnitude of maximum emissions, but the timing is more uncertain, with a range of 18 years for SSP1-2.6 (2015-32) and 15 years for SSP2-4.5 (2043-57); this is due in part to the broader peak and decline of the CO 2 curve in the lower concentration scenarios.For the higher emissions scenarios, SSP5-8.5 and SSP3-7.0particularly, there is a greater spread in the maximum emissions rate; for SSP3-7.0emissions climb right through to 2100, so the year of highest emissions occurs uniformly in the last few years of the simulation, whereas in SSP5-8.5 emissions peak and decline, so there is inevitably uncertainty in the year in which the peak occurs, from 2086 to 2100.

2) CUMULATIVE FOSSIL FUEL EMISSIONS
The cumulative ESM FF emissions are compared with the emission rates from the IAMs (dashed colored lines) in Fig. 9b, with their cumulative totals from 1850 to 2100 in Table 5.The ESMs' cumulative fossil fuel emission rates are depicted individually in Fig. S4.The IAMs' cumulative totals were constructed with the FF emissions data from 1850 to 2014 (Gilfillan et al. 2019).
For the two higher emission scenarios, SSP5-8.5 and SSP3-7.0,although the IAMs' FF emission rates increasingly drop below the ESMs' mean toward the end of the century, periods of higher emissions earlier mean that their cumulative totals are in very close agreement.The IAMs also slightly underestimate the FF rates of the lower emission scenarios compared to the ESMs, but the IAM cumulative total lies within the ESM spread for all scenarios, except SSP4-6.0 for the reasons described above.The ESM mean cumulative emissions for SSP4-3.4match almost exactly those constructed from historical data-based values and the IAM and are almost identical to the those from SSP1-2.6, although all nine ESMs performed the latter while only four performed SSP4-3.4.The atmospheric CO 2 is slightly higher in SSP4-3.4than SSP1-2.6,but the former has much higher crop and pasture area (Fig. 1c), so the impact of land use brings its compatible emissions down in line with those of SSP1-2.6.Generally, where the ESMs disagree with the IAMs, the cumulative emissions from ESMs are higher, although UKESM is the exception; for SSP3-7.0 and SSP5-8.5 UKESM's cumulative emissions are lower than the emissions from IAMs by a small amount.
3) TIMING OF NET-ZERO EMISSIONS REQUIRED FOR 1.58 AND 2.08C OF WARMING SSP1-1.9 has lower radiative forcing than all of the RCPs of CMIP5 and was introduced as a result of the Paris Agreement of 2015 to explore the potential for limiting warming to 1.58C above the preindustrial level (Rogelj et al. 2018).Based on RCP2.6,SSP1-2.6 was designed to be consistent with warming of the order of 2.08C or lower.The assumptions underlying both of these low warming scenarios include extensive investment in sustainability and a reduction in fossil fuel use, leading to net zero fossil fuel emissions during the second half of the twenty-first century, and negative emissions subsequently.Once such an ambitious warming limit has been decided upon, a further decision to make is the target year by which emissions should reduce to zero to achieve this.We will now examine the time frame within which the models predict fossil fuel emissions must drop to zero to adhere to the CO 2 pathway, and the change in global mean temperature the models exhibit under each.
Figure 9a shows that in all ESMs, the emission rate drops to below zero for SSP1-1.9,SSP1-2.6, and SSP4-3.4,the three SSPs with the lowest CO 2 forcing, as well as SSP5-3.4over, the overshoot scenario branching from SSP5-8.5 at 2040 with intensive mitigation.Table 6 shows the year in which the fossil fuel emissions rate first reduced to zero, when smoothed by an 11-yr running mean.
For SSP1-1.9, emissions first become zero in 2056 for the IAM; the five ESMs that ran that scenario all require emissions to become negative, starting in 2056 (UKESM) through to 2071 (CNRM).In the case of SSP1-2.6,all nine ESMs see emissions becoming negative during the period 2076-86, compared to 2076 for the IAM.Only four ESMs ran SSP4-3.4,all requiring a reduction to zero from 2080 to 2091, enclosing the IAM's value of 2084.The overshoot scenario, SSP5-3.4over,sees emissions in the ESMs becoming zero in 2068, as is true of the IAM, through to 2078.
The temperature response of the ESMs to the SSPs is shown as a model mean in Fig. 11, and individually in Fig. S2.The multimodel-mean warming at 2100 relative to the 1850-99 average ranges from 1.58 to 2.28C for the five models that performed SSP1-1.9, with a mean of 1.98C (Table S1).The mean of the nine ESMs' warming at 2100 under SSP1-2.6 is 2.08C, with a range of 1.18C (NorESM2) to 2.88C (CanESM5 and UKESM); five ESMs exhibit warming at 2100 of 2.08C or less under this scenario.

4) LAND CARBON RESPONSE
The ESMs' land carbon response to the scenarios is the net result of a combination of their natural carbon cycle behavior as the CO 2 and climate change, and the extent to which this is perturbed by the prescribed land-use forcing, and the land management practices implemented within the ESM.In the  6.Year at which compatible fossil fuel emissions reduce to zero, where applicable, calculated from the fossil fuel emission rate smoothed with a Savitzky-Golay filter with an 11-yr window.
The multimodel mean land carbon uptake rates (Fig. 12) increase broadly in line with increasing atmospheric CO 2 (Fig. 1a), although the differing extents to which natural land is replaced by crop and pasture across the scenarios (Figs.1c,d) counteract this in some cases.The high CO 2 scenario SSP5-8.5 has the greatest NBP.The SSP3-7.0, however, with second highest CO 2 has the greatest loss of forest that is replaced by the greatest pasture area and the second largest crop area of all SSPs.Additionally, SSP3-7.0 has larger contributions from non-CO 2 greenhouse gases, which do not stimulate land carbon uptake, but reduce terrestrial carbon via positive climate feedback.These factors contribute to the lowering of model mean NBP under SSP3-7.0 to below that of SSP4-6.0 (although only three ESMs performed the latter, compared to all nine for SSP3-7.0).The lowest NBP is simulated in response to SSP4-3.4; it has the third lowest CO 2 pathway, but also features the highest crop area of all scenarios, as well as fourth highest pasture area, which brings its FF emissions down approximately in line with SSP1-2.6 [section 3b(2)].
Although the model mean exhibits some interscenario spread in the time series of NBP, there is very little evidence of this in some models (ACCESS, IPSL, UKESM), while NBP in other models is very scenario dependent (CanESM5, MIROC) (Fig. S4).Figures 13b and 13c resolve the change in land carbon store (Fig. 13a) into change in vegetation carbon, and soil plus litter carbon, respectively.Figures S6-S9 show the individual ESMs' changes in land, vegetation, soil, and litter carbon, respectively.The change in vegetation carbon is responsible, primarily, for the evolution of land carbon throughout the twenty-first century, with model mean change in soil carbon differing by less than 100 GtC across scenarios, as is true for the ESMs individually with the exception of CanESM5 (Fig. S7).
The ESMs' natural carbon cycle behavior can be characterized by the feedback parameters introduced by the original Coupled Climate-Carbon Cycle Model Intercomparison Project (C4MIP) of the mid-2000s (Friedlingstein et al. 2006).The C4MIP project is ongoing and forms one of the model intercomparison projects of CMIP6; it includes a set of idealized experiments designed to allow quantification of carbon-climate and carbon-concentration feedbacks.The idealized 1pctCO2 experiments, in which CO 2 rises at 1% annually to 4 times the preindustrial level, have been performed by fully coupled, biogeochemically coupled, and radiatively coupled variants of 11 CMIP6 ESMs, including all 9 in the present study.This combination of experiments allows the carbon-concentration feedback (denoted b L and b O for land and ocean, respectively, in units of GtC ppm 21 ) and the carbon-climate feedback (g L and g O , GtC 8C 21 ) to be quantified.The feedback parameters help to explain the differences in the models' responses to the SSPs forcings.The greater the b L , the greater the land uptake of carbon per unit increase in atmospheric CO 2 , representing a negative feedback of the carbon cycle-climate system, slowing the accumulation of CO 2 in the atmosphere.A larger, negative g L results in a greater loss of carbon from land per unit warming (a positive feedback, increasing atmospheric CO 2 ).The carbon-climate and carbonconcentration feedbacks of the nine ESMs included in the present study are included in Table A1 of Arora et al. (2020) at 2 times and 4 times preindustrial CO 2 .
Among all models, CanESM5 yields the highest land carbon uptake under SSP5-8.5,leading to higher diagnosed FF emissions than all other ESMs for that scenario.The reason for this is an increase in the strength of CO 2 fertilization effect in CanESM5 compared to its predecessor CanESM2; it has the third highest b L (1.28 GtC ppm 21 at 4 times CO 2 ) of the 11 models in Arora et al. (2020), for which the mean is 0.97 GtC ppm 21 .CanESM5 does not contain an explicit representation of nutrient constraints on land photosynthesis; the strength of the CO 2 fertilization effect over land is determined by a photosynthesis downregulation parameterization.Arora and Scinocca (2016) attempt to constrain the model parameter that determines photosynthesis downregulation, using the historical carbon budget, the amplitude of the annual cycle of globally averaged CO 2 and its trend over the historical period.Yet despite this exercise, the response of land carbon uptake in CanESM5 to increasing CO 2 is larger than other models at high CO 2 .Further, of CNRM has a carbon-concentration feedback comparable to CanESM5 (1.36 GtC ppm 21 at 4 times CO 2 ), contributing to its strong land carbon uptake.It differs from CanESM5 in having a fairly strong carbon-climate feedback (283.11GtC K 21 compared with CanESM5's 115.95 GtC 8C 21 and the multimodel mean of 245.07 GtC 8C 21 at 4 times CO 2 ).As the planet gets progressively warmer under the high CO 2 scenarios, CNRM's positive carbon-climate feedback dampens carbon uptake in a way that does not happen in CanESM5.MIROC also features strong gains in land carbon under SSP5-8.5 and SSP3-7.0,406, and 368 GtC, respectively, consistent with it having the third strongest carbon-concentration feedbacks (1.12 GtC ppm 21 compared with the ESM mean of 0.97).It also has a stronger than average climate-carbon feedback, of 269.57GtC 8C 21 , compared with the ESM mean of 245.07, but MIROC has a lower than average climate sensitivity (Arora et al. 2020) leading to warming in SSP5-8.5 of only around 4.58C.Of the nine ESMs in this study, MIROC, CNRM, and CanESM5 are the three with the largest fraction of emissions taken up by land in Fig. 4b of Arora et al. (2020).
UKESM and IPSL exhibit comparatively weak land carbon uptake, with little dependence on CO 2 pathway due to their relatively weak carbon-climate feedbacks (238.4 and 28.67 GtC 8C 21 , respectively) and carbon-concentration feedbacks of 0.7 and 0.62 GtC ppm 21 .
Since both NorESM2 and CESM2 employ the land carbon cycle model CLM5, their responses in the scenarios are very similar (Fig. S6).They see, respectively, a rise in land carbon of 323 and 305 GtC for SSP5-8.5, and, at the lower end, 189 and 203 GtC for SSP1-2.6.This is consistent with their below average land carbon-concentration feedback.These models have different atmospheric components, and CESM2 sees greater warming than NorESM2, but with much lower than average climate-carbon feedbacks, there is little additional loss of carbon in CESM2 compared to NorESM2, even for SSP5-8.5 (Figs.S6-S9).Arora et al. (2020) discuss some of the differences between models that affect their feedback parameters.The CO 2 fertilization effect, defined as the change in GPP per unit change in atmospheric CO 2 concentration (GtC yr 21 ppm 21 ) contributes to the carbon-concentration feedback; this is higher in the three ESMs in Arora et al. (2020), which include dynamic vegetation, including MPI, and UKESM.However, in both of these models, the turnover time of soil and vegetation carbon are lower than average.CNRM on the other hand, features strong land carbon uptake due to its relatively high turnover times, rather than a particularly strong CO 2 fertilization effect.Another differentiator is the fraction of the change in GPP (relative to the preindustrial) that gets converted to net primary production (NPP); this is a form of carbon use efficiency denoted DCUE to highlight that it is the change in GPP above the preindustrial that is under consideration.CanESM5 has a higher DCUE than all other models, which contributes to its relatively high land carbon uptake.ACCESS sees the lowest gain in land carbon throughout the scenarios because it has the lowest DCUE as well as a weak CO 2 fertilization effect.These combine to give it the lowest carbon-concentration feedback of the nine ESMs in this study.
The idealized experiments used to determine the carbon feedback parameters featured fixed, preindustrial nitrogen deposition and fertilizer.The SSP experiments use timevarying nitrogen deposition, providing an additional driver of land carbon uptake in those ESMs that feature explicit carbonnitrogen interactions, limiting carbon uptake when plantavailable nitrogen is scarce; this is parameterized in others (Tables 1-3 and supplemental section 1).Some models have the ability to diagnose the reduction in land carbon uptake due to the scarcity of nitrogen and other nutrients within a simulation; for example, in UKESM, NPP is calculated before and after N limitation is taken into account, each time step.In the FUN model used in CLM5 within CESM2 and NorESM2, the ''carbon cost'' of nitrogen uptake in units of assimilable carbon is calculated; the higher the carbon cost of nitrogen, the greater the nitrogen limitation (Fisher et al. 2019).However, assessing the full impact of nitrogen limitation on land carbon stores would require spinning up the carbon-nitrogen and carbononly ecosystems in separate simulations, followed by parallel piControl, historical, and scenario experiments, which is prohibitively expensive.
The SSPs' land-use trajectories, forcing the myriad land-use schemes within the ESMs, superimpose another layer of complexity upon the models' biophysical and biological processes, contributing to the changes in land, soil and vegetation carbon stores seen in Figs.S6-S8.The discussion on the ESMs' representation of crop harvest in section 3a(4) is instructive when considering these changes.However, as discussed in section 3a(3), having the no-Lu variants of each SSP scenario would be invaluable in determining the exact role of land use on the land carbon sink by comparing the net change in land carbon storage between the two.

5) OCEAN CARBON RESPONSE
There is much more agreement between ESMs in their ocean carbon response to the scenarios than is the case for land carbon, as shown in Figs.14a and 14b, depicting the air-to-sea flux and change in ocean carbon store, respectively; the same are shown for the models individually in Figs.S10 and S11, respectively.The air-to-sea flux averaged across a portion of the piControl was subtracted from that of the scenarios for all ESMs as in some cases the piControl ocean was not fully in equilibrium, with a small, but not negligible, flux to or from the atmosphere.The separation between scenarios is significant in all ESMs, and their cumulative uptake totals are in close agreement (Table S5).ACCESS gains slightly more carbon than all other ESMs over the historical period, which continues throughout the SSP scenarios leading to the greatest uptake at 2100, relative to 1850, of all ESMs for its four scenarios.This is consistent with it having the largest carbon-concentration feedback (0.9 GtC ppm 21 ) of all the ESMs included in this study for which it is known [the other eight range from 0.70 to 0.78 GtC ppm 21 ; Table A1 in Arora et al. (2020)].However, its negative climate-carbon feedback (223.75GtC 8C 21 ) is also largest in magnitude of all ESMs, the others ranging from 29.38 to 222.25 GtC 8C 21 , promoting greater carbon loss per 18C of warming than other ESMs, offsetting partly the gain in carbon.
All models exhibit either a stabilization or a decline in ocean carbon uptake toward the end of the twenty-first century in all scenarios except for SSP3-7.0,where carbon uptake continues to increase.The reason for this is that the transient trend in ocean carbon uptake strongly depends on the trend of atmospheric CO 2 or dCO 2 /dt (Tjiputra et al. 2014).The CO 2 trend in SSP3-7.0 continues to increase by 2100, whereas it stabilizes or declines in other SSPs.
The ESMs' ocean carbon response to the scenarios is therefore much less diverse than is true for land carbon, as was the case CMIP5 (Jones et al. 2013).Although there is diversity in the nature and complexity of ecosystems included in the ocean carbon cycle models, and to an extent the same is true of the carbonate chemistry and air-to-sea flux parameterizations, they are less conflicted than the terrestrial carbon components, which need to incorporate the effects of human interaction with land.

6) AIRBORNE FRACTION
The cumulative airborne fraction of fossil fuel emissions for the SSPs are shown in Fig. 15 and individually in Fig. S12.
From the initial large model differences at the end of the historical simulation (Fig. 7), there is quite strong agreement among ESMs that by the second half of the twenty-first century the airborne fraction will increase for the higher emission scenarios, SSP5-8.5 and SSP3-7.0, and reduce for the lower emission SSP1-1.9 and SSP1-2.6.Indeed, for the latter two, the models are in extremely close agreement, with the model mean CAF reducing by 0.11 over the last half of the century in both scenarios.In the higher emissions scenarios, the sinks cannot keep pace with the amount of CO 2 being released, and the fraction left in the atmosphere increases.For the lower emission scenarios, the emission rate is lower than the combined sink strength, and therefore the fraction of emissions left in the atmosphere reduces.
For the intermediate scenarios, SSP4-3.4,SSP4-6.0, and SSP5-3.4over, the CAF tends to reduce, but at a lower rate than is the case for SSP1-1.9 and SSP1-2.6, with some individual models being almost constant over the last few decades.The overshoot scenario, SSP5-3.4over,sees a slightly greater decline in CAF than SSP1-1.9 and SSP1-2.6, at 0.14 over the last 50 years.The intermodel spread is such that the highest CAF is typically around 0.15 to 0.20 greater than the lowest CAF at 2100.CNRM has the lowest CAFs almost for the duration of all scenarios.

Conclusions
For nine state-of-the-art, comprehensive Earth system models, we present the fossil fuel CO 2 emission rates compatible with the CO 2 concentration pathway used to drive them over the historical period, and under the eight Shared Socioeconomic Pathway scenarios out to 2100.
The model-mean time series of simulated historical fossil fuel CO 2 emissions compares very well with data-based estimates; the mean cumulative total from 1850 to 2014 (394 6 59 GtC) is extremely close to the estimate from the global carbon budget 2019 (400 6 20 GtC); only two models, however, lie within the observed estimate's uncertainty range.With the exception of CNRM and ACCESS, the models that are outside of the observed uncertainty range underestimate cumulative historical emissions.A major contributing factor to this is the stabilization of observed CO 2 concentration in the 1940s that led to lower diagnosed FF emissions in the ESMs compared to reality, which saw emissions continue to rise.
Good agreement between the observed rate of historical fossil emissions and those from ESMs provides confidence in their ability to simulate the net carbon uptake by land and ocean reliably over the historical period.Good agreement between fossil emissions diagnosed from the immensely complex, process-based ESMs and the emissions calculated by the much simpler IAMs for the future scenarios provides confidence in our understanding of the coupled carbon cycle and climate system.Indeed, agreement between the two for the SSP experiments is quite strong, when looking at both time series and the cumulative totals of compatible fossil fuel CO 2 emissions.If this were not the case it would undermine credibility in the sequence of influence from the scenarios of population change and socioeconomic variables driving the IAM, the resulting CO 2 concentration pathway, and the physical and biogeochemical changes simulated by the ESM in response to it, as well as associated climate impacts.However, they are sufficiently close that we can have confidence that the CO 2 scenarios forcing the complex models are a realistic expression of the evolving worlds described by the Shared Socioeconomic Pathway scenarios.The fact that the IAM emissions are slightly lower than the ESMs' average would imply that the land and ocean uptake in those models is perhaps lower than that of the ESMs.
The compatible emissions time series presented here provide an uncertainty range as well as context to the single CO 2 emissions rate generated by the IAM for each SSP.Further, with the year by which society should aspire to reduce global emissions to zero gaining ever more attention, the ESMs provide an uncertainty range around this event.
As with CMIP5, the contrast between the land and the ocean carbon responses to future forcing is marked; both generations of models exhibit a strong lack of agreement in the response of the land carbon sink to the diverse range of CO 2 scenarios, with much greater agreement for ocean carbon.One primary reason land models diverge more than ocean models in terms of their future projections of carbon uptake is that ocean carbon uptake is determined almost entirely by well understood physical and chemical processes; over land, carbon uptake is determined exclusively by biological processes that are not as well understood.In addition, uncertainty over land is contributed to by several other processes including land-use change and land management practices (e.g., treatment of changing crop and pasture area, fertilizer application, irrigation, harvesting of biomass), as well as differences in nutrient cycling (nitrogen, phosphorus), all of which affect the land carbon sink.
ESMs are ever-evolving, with some already including the land carbon impact of fire-vegetation interactions.The release of CO 2 and methane from thawing permafrost (as already occurs in CESM2) as well as methane from wetlands will likely be included in future ESMs, and the development of interactive methane emissions with a 3D atmospheric CH 4 tracer, in parallel with the CO 2 equivalent, offers further opportunities to investigate feedbacks between the carbon cycle, the climate and the wider Earth system.A continuing goal in the coming decade should be to better constrain the land carbon response, despite these myriad differences.
The ''E2'' paired historical and hist-noLu experiments provide the best estimate of land-use and land-cover change emissions, E LUC , in the historical experiment.The seven ESMs that have provided output from both experiments to the CMIP6 archive all underestimate land-use emissions compared to the Unauthenticated | Downloaded 05/20/24 03:33 AM UTC two reference datasets from the global carbon budget 2019 used in Fig. 3, to a greater or lesser extent, despite the GCB2019 data not accounting for the loss of additional sink capacity included in the ESMs.The ESMs at the edges of the model spread in E FF and E LUC (Figs. 3b,d) suggest an inverse relationship between the two; greater land-use emissions are paired with fewer fossil fuel emissions.This relationship suggests that in the absence of the -noLu variant of the SSP scenarios, the magnitude of an ESM's historical LULCC emissions from paired experiments could potentially be used, with caution, to provide a way of qualitatively interpreting its FF emissions.If, for example, an ESM yields very high (low) historical paired-experiment LULCC emissions, its fossil fuel emissions diagnosed from future scenarios could be considered excessively low (high).
Using the E2 experimental design to allow comparison of FF, LULCC and FF 1 LULCC emissions between models for future scenarios, as well as for the historical period, would allow us to examine how future land use, accounting for loss of additional sink capacity and regrowth on abandoned agricultural land, impact FF emissions.It would also allow comparison between models of the response of the natural carbon cycle to realistic, rather than idealized scenarios, due primarily to changing environmental conditions rather than to the direct effects of land use on land carbon, as discussed in section 3a(3) in relation to Fig. 3f.
Aerosols, well-mixed greenhouse gases and ozone all influence the evolving climate and, therefore, land and ocean carbon uptake, and fossil emissions.The role of non-CO 2 forcings on the carbon cycle, future carbon budgets, and on the spread of the land carbon response among ESMs in particular, is an area that would benefit from further attention in the coming years.As models incorporate new processes such as interactive diffuse:direct fraction of photosynthetically active radiation, online calculation of atmospheric methane concentration from emissions, interactive nitrogen deposition as a product of atmospheric chemistry processes, as well as downgrading of carbon uptake by plants due to ozone damage, this will likely become increasingly important.
ESMs can be driven with prescribed atmospheric CO 2 concentration (C driven), as in the experiments described here, or provided with time varying fossil fuel CO 2 emissions (E driven) from which the model calculates the atmospheric CO 2 concentration.In C-driven mode, the diagnosed fossil fuel flux is the residual of the top-level carbon fluxes between the atmosphere and the land and the ocean, so the uncertainty associated with every physical, biological, and biogeochemical process that contributes to the partitioning of carbon between land and ocean accumulates in the residual FF emissions.Conversely, in E-driven mode with prescribed fossil fuel CO 2 emissions, the same uncertainties that affect the partitioning between land and ocean feed into the changing atmospheric CO 2 burden.The diagnosed fossil fuel emissions, or the changing atmospheric CO 2 burden, therefore, tend to act as integrators of model spread when comparing multiple ESMs.A useful complement to the present paper would be a similar study examining the emissions-driven experiments of the CMIP6 Earth system models.

FIG. 1 .
FIG. 1. Scenarios of (a) atmospheric CO 2 concentration (parts per million), (b) global population (billions of people), and area (millions of hectares) of the land surface covered by (c) crop and pasture and (d) forest, of the eight SSPs.Historical CO 2 concentration is shown in the inset in (a).

FIG. 2
FIG. 2. (a) Time series of compatible fossil fuel emissions (GtC yr 21 ) from 1850 to 2014 simulated by the nine ESMs.Solid black line is the model mean.The dashed black and white line is the data-based fossil fuel rate from the Gilfillian et al. (2019) as cited in the GCB2019 (Friedlingstein et al. 2019), with decadal mean emission rates from the 1960s to the 2010s and associated uncertainty, indicated by the white error bars, from GCB2019.(b) Cumulative fossil fuel emissions (GtC) for the same period.
FIG. 3. (left) Time series and (right) cumulative totals of CO 2 emissions from seven ESMs with observational fossil fuel emissions from Gilfillan et al. (2019) and land-use emissions from the GCB2019, which also provides the uncertainty range in cumulative totals at 2014 (Friedlingstein et al. 2019).(a) Compatible fossil fuel emissions (GtC yr 21 ) calculated for the historical experiment.(b) Cumulative fossil fuel emissions from the historical experiment.(c) Landuse emissions calculated as land carbon uptake (cumulative NBP) in hist-noLu minus that in historical.(d) Cumulative land-use emissions.(e) Combined fossil fuel and land-use emissions.(f) Cumulative total of combined fossil fuel and land-use emissions.Note that (c) to (f) contain two data-based estimates of E LUC .The thick black curve is the dataset of land-use emissions provided under the C4MIP protocols for driving the emissions-driven historical esm-hist experiment for models that do not have the ability to simulate LUC emissions interactively with the atmosphere (Jones 2016).These are derived from Hansis et al. (2015) and Houghton et al. (2012).The thin black curve is E LUC from Friedlingstein et al. (2019), the cumulative total calculated using Table8therein.

FIG. 4 .
FIG. 4. Time series of net biosphere production (GtC yr 21 ) from 1850 to 2014 calculated by the ESMs.Data-based estimates of decadal mean net land carbon uptake for the1960s, 1970s, 1980s,  1990s, and 2000s and from 2009  to 2018 are overlaid as squares with uncertainty range, from Table6ofFriedlingstein et al (2019).

FIG. 5
FIG. 5. (a) Initial land carbon stores in all ESMs at 1850, partitioned into vegetation, litter, and soil carbon.(b) Change in land carbon from 1850 to 2014 calculated as the cumulative nbp diagnostic.(c) Change in vegetation carbon and (d) change in 1litter carbon over the same period.The values in (c) and (d) are calculated from the cVeg and cSoil 1 cLitter diagnostics, respectively.
FIG. 7. (a) Time series of air-to-sea flux of CO 2 from the nine ESM's historical experiments.Data-based estimates of decadal mean fluxes for the 1960s, 1970s, 1980s, 1990s, and 2000s and from 2009 to 2018 are overlaid as squares with an uncertainty range.(b) Change in ocean carbon store calculated as the cumulative total of the air-to-sea flux shown in (a).The square at 2014 indicates the data-based estimate of cumulative carbon uptake from 1850 to 2014.All observations are from Table 6 of Friedlingstein et al. (2019).

FIG. 6
FIG. 6.(a) Harvest carbon flux from the historical simulation of five ESMs.The data are annual means calculated from one ensemble member, with the average fHarvest from at least a 10-yr section of the piControl subtracted.(b) Cumulative totals of harvest fluxes shown in (a).

FIG. 8 .
FIG. 8. Time series of airborne fraction of fossil fuel emissions from the 9 ESM's historical simulation.Black squares and uncertainty ranges are observations-based estimates from Table 6.1 of Ciais et al. (2013).

FIG. 9
FIG. 9. (a) Time series of compatible fossil fuel emissions (GtC yr 21 ) for the eight SSP scenarios simulated by the nine ESMs.Solid lines are the model mean, the dashed lines are the fossil fuel rate generated by the IAMs that provided the CO 2 concentration pathways used to drive the ESMs.Shaded areas indicate model spread for each scenario.(b) Cumulative totals of fossil fuel emissions (GtC) shown in (a).

FIG. 10 .
FIG. 10.Maximum rate of FF emissions (GtC yr 21 ) experienced during each SSP by each ESM vs year in which this occurred.Symbols denote SSP and ESMs are distinguished by color, with IAM in black.

FIG. 13 .
FIG. 13.Net change (GtC) in (a) land carbon store, (b) vegetation carbon, and (c) soil 1 litter carbon for the SSPs.Thick lines are the multimodel means; individual models are indicated by the thin lines.Calculated from the cumulative nbp, and change in cVeg, and change in cSoil 1 cLitter diagnostics, respectively.
FIG. 14.(a) Air-to-sea flux of CO 2 for the eight SSPs.Multimodel mean for each SSP is shown as the solid line, with model spread indicated by the shaded region around it.(b) Change in ocean carbon store for the SSPs.Thick lines are model means; thin lines are individual models.

FIG. 15 .
FIG. 15.Cumulative airborne fraction of emitted fossil fuel CO 2 from 2015 to 2100 for the eight SSPs simulated by the ESMs.Colored lines depict individual models, with the multimodel mean as the black line.

TABLE 3 .
Summary of MPI-ESM1.2-LR,NorESM2-LM, and UKESM1-0-LL.See supplementary information for detailed model descriptions and explanation of acronyms used in the table.